clear all;
clc;

% ----------------------- Analysis Parameters -----------------------------

pol = -1;
noaxons = 200;

patientgeom = load('femspine_patientgeom.txt');
cordgeom = load('patient_cordgeom.txt');

epidur  = 1; % 0 = intradural, 1 = epidural
subloc  = 1; % 1 = 1 mm above cord; 2 = 1 mm below dura
th_elec = 0; % 0, -10, -20
thtag = ['Lateral Deviation = ',num2str(th_elec),'^o'];

patients = 1:5;
nopat = length(patients);

dc_th   = zeros(noaxons,nopat);
dr_th   = zeros(noaxons,nopat);
ra_elec = zeros(1,nopat);

% experimental data
if(epidur==0)
    isens_exp = [0.20,0.35,0.70,0.55,0.80];
    ipain_exp = [0.50,0.50,1.50,1.40,3.10];
else
    isens_exp = [1.30,19.0,7.00,5.00,9.00];
    ipain_exp = [2.20,24.5,19.0,5.00,19.0];    
end

% ----------------------- Importing Data ----------------------------------

dcderm_data = cell(1,nopat);

for ii = 1:nopat
    
    patnum = patients(ii);

    % location of electrode
    if(epidur==1)
        eletag = 'epi_';
        fdist_elec = cordgeom(patnum,8);
    else
        eletag = 'sub_';
        if(subloc==1)
            fdist_elec = cordgeom(patnum,9);
        elseif(subloc==2)
            fdist_elec = 1-cordgeom(patnum,9);
        else
            return;
        end
    end
    loctag = ['p',num2str(patnum),'_d',strrep(num2str(fdist_elec),'.','p'),...
              '_t',strrep(num2str(th_elec),'-','n')];

    p_tag = num2str(patnum);

    prefix = ['neg_',eletag,'bp_'];
    suffix = [loctag,'.txt'];

    cd(['patient',p_tag,'_thresholds']);
    dcthresh = load([prefix,'dc_act_',suffix]);
    drthresh = load([prefix,'dr_act_',suffix]);
    cd ..;

    cd(['patient',num2str(patnum),'_potentials']);
    iunit = load([eletag,'bp_current_',loctag,'.txt']);
    rtiss = 1/iunit/1e3;
    cd ..;
    
    dc_th(:,ii) = dcthresh(:,1);
    dr_th(:,ii) = drthresh(:,1);
    ra_elec(ii) = rtiss;
    
    dcth_derm = zeros(10,10);
    dcth_derm(:) = -dcthresh(1:(noaxons/2),1);
    dcderm_data{ii} = dcth_derm/rtiss;
    
end

aa = 1:10;
bb = ( aa'*ones(1,10) )';

RA = ( ra_elec'*ones(1,noaxons) )';
dc_th = dc_th./RA;
dr_th = dr_th./RA;

% ----------------------- Processing Data ---------------------------------

ratio_exp = ipain_exp./isens_exp;

thtag = ['Lateral Deviation = ',num2str(th_elec),'^o'];
t8_vertlv  = {'T9','T10/11','T12','L1','L2/3','L3/4','L5','S1/2','S2/3','S4/5'};

noy = 10;
for ii = 1:2
    
    xx = bb(:);
    yy = dcderm_data{ii}(:);
    figure;
    
    mindcth = min(yy);
    prepain = ratio_exp(ii)*mindcth;
    
    subplot(noy,1,1:(noy-1));
    hold on;
    plot(xx,yy,'k.','MarkerSize',30);
%     plot(aa,mindcth*ones(size(aa)),'k--','LineWidth',3);
%     plot(aa,prepain*ones(size(aa)),'k-','LineWidth',3);
    ttl = ['Patient ',num2str(ii),', ',thtag];
    title(ttl,'FontSize',30);
    ylabel('Stimulation Amplitude (mA)','FontSize',30);
    hold off;
    xlim([0,11]);
    ylim([0,20]);
    L1 = 'DC Stimulation Threshold';
    L2 = 'Predicted Sensory Threshold';
    L3 = 'Predicted Pain Threshold';
    legend(L1,'Location','NE');
    set(gca,'XTick',[],'FontSize',26);
    set(gca,'YScale','Log');
    subplot(noy,1,noy);
    xlim([0,11]);
    ylim([0,1]);
    for jj = 1:length(t8_vertlv)
        nolet = length(t8_vertlv{jj});
        text(jj-0.175*nolet/2,0.5,t8_vertlv{jj},'FontSize',30);
    end    
    axis off;
    
end
